### regression function

regs_int <- function(panel_only = TRUE, d) {
  
  cols = c("dv", "year", "treat", "mod",
             "kecid", "pidlink")

  d = d[complete.cases(d[,cols]),]
  
  if (panel_only) {
    
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
    
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)    
  }

  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  ### clustered standard errors
  d$kecid = factor(d$kecid)
  
  m = lme4::lmer(dv ~ year*treat*mod + (1|kecid), 
         data=d)
  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 4)
  df$SE = round(df$SE, 4)
  df$t = round(df$t, 3)
  df$p = round(df$p, 3)

  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d))
}


regs_int_cov <- function(panel_only = TRUE, shalat = TRUE, how_religious = TRUE, d) {
  
  if (panel_only) {
    
    cols = c("dv", "mod", "year", "treat", 
             "female", "age", "edu", "how_religious", "shalat3pt", "work_lastyear", "marital", 
             "javaisland", "census10_median_age", "census10_prop_male", "census10_prop_working", 
             "census10_median_edu",
             "census10_prop_muslim", "census10_diversity_lang", 
             "elec14_kec_jokowi_share", "elec19_islamist_share", "elec19_jokowi_share",
             "n_mo_07", 
             "nmprov_mosque", "nmkab_mosque", "nmkec_mosque", "kecid", "pidlink")
    
    d = d[complete.cases(d[,cols]),]
    
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
    
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)
    
  }
  
  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  ### clustered models
  d$kecid = factor(d$kecid)
  
  if (shalat & how_religious) { 
    m = lme4::lmer(dv ~ year*treat*mod +
               female + age + edu + how_religious + shalatlab3pt + work_lastyear + marital + 
               javaisland + census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
               census10_prop_muslim + census10_diversity_lang + 
               elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
               ln_n_mo_07 + (1|kecid), 
             data=d)
  } 
  
  if (!shalat) {
    m = lme4::lmer(dv ~ year*treat*mod + 
               female + age + edu + how_religious + work_lastyear + marital + 
               javaisland + census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
              census10_prop_muslim + census10_diversity_lang + 
              elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
              ln_n_mo_07 + (1|kecid), 
           data=d)
  }

  if (!how_religious) {
    m = lme4::lmer(dv ~ year*treat*mod + 
             female + age + edu + shalatlab3pt + work_lastyear + marital + 
             javaisland + census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
             census10_prop_muslim + census10_diversity_lang + 
             elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
             ln_n_mo_07 + (1|kecid), 
            data=d)
  }
  
  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 5)
  df$SE = round(df$SE, 5)
  df$t = round(df$t, 3)
  df$p = round(df$p, 3)

  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d))
}



regs_basic <- function(panel_only, d) {
  
  cols = c("dv", "year", "treat")
  
  d = d[complete.cases(d[,cols]),]
  
  if (panel_only) {
  
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
  
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)
  }
  
  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  
  ### clustered standard errors
  d$kecid = factor(d$kecid)
  
  m = lme4::lmer(dv ~ year:treat + year + treat + (1|kecid), 
                    data=d)
  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 4)
  df$SE = round(df$SE, 4)
  df$t = round(df$t, 3)
  # df$p = round(df$p, 3)
  
  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d, m_cl = m_cl))
}


regs <- function(panel_only, shalat = TRUE, how_religious = TRUE, d) {
  
  if (panel_only) {
    
    cols = c("dv", "year", "treat", 
             "female", "age", "edu", "how_religious", "shalat3pt", "work_lastyear",
           "marital", "javaisland", 
           "census10_median_age", "census10_prop_male", "census10_prop_working", "census10_median_edu",
           "census10_prop_muslim", "census10_diversity_lang", 
           "elec14_kec_jokowi_share", "elec19_islamist_share", "elec19_jokowi_share",
           "n_mo_07", 
           "nmprov_mosque", "nmkab_mosque", "nmkec_mosque", "kecid", "pidlink")

    d = d[complete.cases(d[,cols]),]
  
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
  
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)

  }

  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  
  ### basic models
  d$kecid = factor(d$kecid)
  
  if (shalat & how_religious) {

      m = lme4::lmer(dv ~ year:treat + year + treat + 
               female + age + edu + how_religious + shalatlab3pt + work_lastyear + marital + 
               javaisland + 
               census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
               census10_prop_muslim + census10_diversity_lang + 
               elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
               ln_n_mo_07 + (1|kecid), 
             data=d)
  } 

  if (!shalat) {
    m = lme4::lmer(dv ~ year:treat + year + treat + 
                     female + age + edu + how_religious + work_lastyear + marital + 
                     javaisland + 
                     census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
                     census10_prop_muslim + census10_diversity_lang + 
                     elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
                     ln_n_mo_07 + (1|kecid), 
                   data=d)
  }
  

  if (!how_religious) {
    m = lme4::lmer(dv ~ year:treat + year + treat + 
                     female + age + edu + shalatlab3pt + work_lastyear + marital + 
                     javaisland + 
                     census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
                     census10_prop_muslim + census10_diversity_lang + 
                     elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
                     ln_n_mo_07 + (1|kecid), 
                   data=d)
  }
  
  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 4)
  df$SE = round(df$SE, 4)
  df$t = round(df$t, 3)
  df$p = round(df$p, 3)

  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d))

}


regs_matched <- function(panel_only = FALSE, d = dat) {
  
  if (panel_only) {
    
    cols = c("dv", "year", "treat", "distance", "kecid", "pidlink")
    
    d = d[complete.cases(d[,cols]),]
    
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
    
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)
    
  }
  
  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  
  ### clustered se models
  d$kecid = factor(d$kecid)
  
  m = lme4::lmer(dv ~ year:treat + year + treat + distance + year:distance + (1|kecid),
           data=d)

  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 4)
  df$SE = round(df$SE, 4)
  df$t = round(df$t, 3)
  df$p = round(df$p, 3)

  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d))
  
}


### plotting function

figs <- function(out, ylab, title, ymax = 100) {
  
  ### predict with clustered SE
  # https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html
  
  pred <- ggpredict(out$m, terms = c("year","treat"))
  pred <- ggpredict(out$m, terms = c("year","treat"),
                    vcov.fun = "vcovCR", vcov.type = "CR0",
                    vcov.args = list(cluster = out$data$kecid))
  
  
  pred$predicted = round(pred$predicted*100,2)
  pred$std.error = round(pred$std.error*100,2)
  p <- ggplot(pred, aes(x = x, y = predicted, fill = group)) +
    geom_bar(stat = "identity",position = position_dodge(.9)) +
    geom_errorbar(aes(ymin=predicted-qnorm(.975)*std.error, 
                      ymax=predicted+qnorm(.975)*std.error), width=.2,
                  position=position_dodge(.9)) + 
    geom_text(aes(label=predicted, y = -3), color="black",
              position = position_dodge(0.9), size = 3) +
    scale_y_continuous(name = ylab,
                       limits = c(-5, 100),
                       breaks = seq(0,100,10)) +
    scale_x_discrete(name = "") +
    theme_bw() +
    theme(legend.position = "bottom",
          axis.text.x = element_text(size = 10)) +
    labs(subtitle = title,
         fill = NULL)
  
  
  print(p)
  
}


### regression coefficient names

change_coefs <- function (vec){
  
  vec[vec == "year2014"] = "Year(2014)"
  vec[vec == "treatNew Mosque(s)"] = "Treatment: Mosques"
  
  vec[vec == "female"] = "Female"
  vec[vec == "age"] = "Age"
  vec[vec == "edu"] = "Education"
  vec[vec == "how_religious"] = "Religiosity"
  vec[vec == "shalatlab3ptless than five"] = "Salat: Less than Five"
  vec[vec == "shalatlab3ptmore than five"] = "Salat: More than Five"
  vec[vec == "work_lastyear"] = "Worked Last Year"
  
  vec[vec == "marital"] = "Married"
  vec[vec == "javaisland"] = "Java Island"
  vec[vec == "census10_median_age"] = "Census: Median Age"
  vec[vec == "census10_prop_male"] = "Census: Prop. Males"
  vec[vec == "census10_prop_working"] = "Census: Prop. Working"
  vec[vec == "census10_median_edu"] = "Census: Median Education"
  vec[vec == "census10_prop_muslim"] = "Census: Prop. Muslims"
  vec[vec == "census10_diversity_lang"] = "Census: Language"
  vec[vec == "elec14_kec_jokowi_share"] = "Jokowi Share (2014)"
  vec[vec == "elec19_islamist_share"] = "Islamist Share (2019)"
  vec[vec == "elec19_jokowi_share"] = "Jokowi Share (2019)"
  vec[vec == "n_mo_07"] = "N Mosques (2007)"
  vec[vec == "ln_n_mo_07"] = "Log(N Mosques in 2007)"
  vec[vec == "sharia_pre08"] = "Existing Sharia Ordinances"
  vec[vec == "sharia_0813"] = "New Sharia Ordinances"
  
  vec[vec == "year2014:treatNew Mosque(s)"] = "Year X Treat"
  
  
  vec[vec == "treatOne New Mosque"] = "Treat: One New Mosque"
  vec[vec == "year2014:treatOne New Mosque"] = "Year X Treat"
  
  vec[vec == "distance"] = "Propensity"
  vec[vec == "year2014:distance"] = "Year X Propensity"
  
  vec[vec == "treatNew Mushalla(s)"] = "Treatment: Mushallas"
  vec[vec == "year2014:treatNew Mushalla(s)"] = "Year X Treat"
  
  vec[vec == "mod"] = "Moderator"
  vec[vec == "year2014:mod"] = "Year X Moderator"
  vec[vec == "treatNew Mosque(s):mod"] = "Treat X Moderator"
  vec[vec == "year2014:treatNew Mosque(s):mod"] = "Year X Treat X Moderator"
  
  vec[vec == "modNU/Muhammadiyah"] = "Moderator (NU/Muhammadiyah)"
  vec[vec == "year2014:modNU/Muhammadiyah"] = "Year X Moderator"
  vec[vec == "treatNew Mosque(s):modNU/Muhammadiyah"] = "Treat X Moderator"
  vec[vec == "year2014:treatNew Mosque(s):modNU/Muhammadiyah"] = "Year X Treat X Moderator"
  
  vec[vec == "modNU"] = "Moderator (NU)"
  vec[vec == "modMuhammadiyah"] = "Moderator (Muhammadiyah)"
  vec[vec == "year2014:modNU"] = "Year X NU"
  vec[vec == "year2014:modMuhammadiyah"] = "Year X Muhammadiyah"
  vec[vec == "treatNew Mosque(s):modNU"] = "Treat X NU"
  vec[vec == "treatNew Mosque(s):modMuhammadiyah"] = "Treat X Muhammadiyah"
  vec[vec == "year2014:treatNew Mosque(s):modNU"] = "Year X Treat X NU"
  vec[vec == "year2014:treatNew Mosque(s):modMuhammadiyah"] = "Year X Treat X Muhammadiyah"

  vec[vec == "shalat3pt_1"] = "Salat: Five Times"
  vec[vec == "shalat3pt_2"] = "Salat: Less than Five"
  vec[vec == "shalat3pt_3"] = "Salat: More than Five"
  vec[vec == "live_village2pt"] = "Reject: Village"
  vec[vec == "live_neighborhood2pt"] = "Reject: Neighborhood"
  vec[vec == "live_room2pt"] = "Reject: House"
  vec[vec == "worship_house2pt"] = "Reject House of Worship"
  vec[vec == "intermarriage2pt"] = "Reject Intermarriage"
  vec[vec == "similar_religion"] = "Prefer Co-Religionist"
  vec[vec == "similar_ethnicity"] = "Prefer Co-Ethnic"
  vec[vec == "trust_religion2pt"] = "Trust Fellow Muslims"
  vec[vec == "trust_ethnicity2pt"] = "Trust Co-Ethnic"
  vec[vec == "helping2pt"] = "Helping Neighbors"
  
  
  return(vec)
}



plot_parallel <- function(start, file) {

  ##### trends for control and treatment
  
  trends = kec %>%
    group_by(year, Treatment) %>%
    dplyr::summarise(meanviol_all = mean(nviol_all, na.rm=TRUE),
                     meanviol_rel = mean(nviol_rel, na.rm=TRUE),
                     meanviol_intrarel = mean(nviol_intrarel, na.rm=TRUE),
                     meanviol_interrel = mean(nviol_interrel, na.rm=TRUE),
                     meanviol_vice = mean(nviol_vice, na.rm=TRUE),
                     prviol_intrarel = mean(prviol_intrarel, na.rm=TRUE),
                     prviol_interrel = mean(prviol_interrel, na.rm=TRUE),
                     prviol_rel = mean(prviol_rel, na.rm=TRUE),
                     prviol_vice = mean(prviol_vice, na.rm=TRUE))
  
  trends = trends[complete.cases(trends[,c("year", "meanviol_all", "Treatment")]),]
  
  p_all = ggplot(data = trends, aes(x = year, y = meanviol_all, group = Treatment)) +
    geom_line(aes(color = Treatment, lty = Treatment), lwd = 1.5) +
    geom_vline(xintercept = 2007, lty = 2) +
    geom_vline(xintercept = 2014, lty = 2) +
    geom_label(aes(x = 2007, y = 50), label = "IFLS-4") + 
    geom_label(aes(x = 2014, y = 50), label = "IFLS-5") + 
    geom_label(aes(x = start+(2013-start)/2, y = 90), label = "Treatment Period\n(Endpoints Included)") +
    annotate(geom = "rect", xmin = start, xmax = 2013, ymin = 0, ymax = 100,
             fill = "gray", alpha = 0.3) +
    scale_x_continuous(name = "Year", 
                       breaks = seq(0, max(trends$year), 2),
                       limits = c(1998,2015)) + 
    scale_y_continuous(name = "Average Number of Incidents") +
    labs(subtitle = "All Incidents",
         color = NULL,
         lty = NULL) + 
    theme_bw() +
    theme(legend.position = "bottom")
  
  p_relprop = ggplot(data = trends, aes(x = year, y = prviol_rel*100, group = Treatment)) +
    geom_line(aes(color = Treatment, lty = Treatment), lwd = 1.5) +
    geom_vline(xintercept = 2007, lty = 2) +
    geom_vline(xintercept = 2014, lty = 2) +
    geom_label(aes(x = 2007, y = 3), label = "IFLS-4") + 
    geom_label(aes(x = 2014, y = 3), label = "IFLS-5") + 
    geom_label(aes(x = start+(2013-start)/2, y = 4), label = "Treatment Period\n(Endpoints Included)") +
    annotate(geom = "rect", xmin = start, xmax = 2013, ymin = 0, ymax = 5,
             fill = "gray", alpha = 0.3) +
    scale_x_continuous(name = "Year", 
                       breaks = seq(0, max(trends$year), 2),
                       limits = c(1998,2015)) + 
    scale_y_continuous(name = "Percentage",
                       limits = c(0, 5)) +
    labs(subtitle = "Percentage of Religion-Related Incidents Relative to Total Incidents",
         color = NULL,
         lty = NULL) + 
    theme_bw() +
    theme(legend.position = "bottom")


  p_relmean = ggplot(data = trends, aes(x = year, y = meanviol_rel, group = Treatment)) +
    geom_line(aes(color = Treatment, lty = Treatment), lwd = 1.5) +
    geom_vline(xintercept = 2007, lty = 2) +
    geom_vline(xintercept = 2014, lty = 2) +
    geom_label(aes(x = 2007, y = .3), label = "IFLS-4") + 
    geom_label(aes(x = 2014, y = .3), label = "IFLS-5") + 
    geom_label(aes(x = start+(2013-start)/2, y = .4), label = "Treatment Period\n(Endpoints Included)") +
    annotate(geom = "rect", xmin = start, xmax = 2013, ymin = 0, ymax = .44,
             fill = "gray", alpha = 0.3) +
    scale_x_continuous(name = "Year", 
                       breaks = seq(0, max(trends$year), 2),
                       limits = c(1998,2015)) + 
    scale_y_continuous(name = "Average",
                       limits = c(0, .45)) +
    labs(subtitle = "Average Number of Religion-Related Incidents",
         color = NULL,
         lty = NULL) + 
    theme_bw() +
    theme(legend.position = "bottom")
  
  
  #### treatment effects
  
  coefs = read.csv(file, header=TRUE)
  
  coefs$low99 = coefs$coef - qnorm(.995) * coefs$se
  coefs$hi99 = coefs$coef + qnorm(.995) * coefs$se
  
  coefs$low95 = coefs$coef - qnorm(.975) * coefs$se
  coefs$hi95 = coefs$coef + qnorm(.975) * coefs$se
  
  coefs$low90 = coefs$coef - qnorm(.95) * coefs$se
  coefs$hi90 = coefs$coef + qnorm(.95) * coefs$se
  
  dvs = c("Object: Live in Same Village", 
          "Object: Live in Same Neighborhood", 
          "Object: Live in Same House",
          "Object: Worship House", 
          "Object: Intermarriage", 
          "Prefer Muslim Candidate", 
          "Prefer Co-Ethnic Candidate",
          "Trust Fellow Muslims", 
          "Trust Co-Ethnics", 
          "Helping Neighbors")
  coefs$dv = dvs
  
  coefs$dv = factor(coefs$dv,
                    c("Helping Neighbors", 
                      "Trust Co-Ethnics", 
                      "Trust Fellow Muslims", 
                      "Prefer Co-Ethnic Candidate",
                      "Prefer Muslim Candidate", 
                      "Object: Intermarriage",
                      "Object: Worship House",
                      "Object: Live in Same House",
                      "Object: Live in Same Neighborhood", 
                      "Object: Live in Same Village"
                    )
  )
  
  coefs$category = c(rep("Outgroup\nRejection", 5),
                     rep("Political\nPreferences", 2),
                     rep("Ingroup\nAttitudes", 3))
  
  coefs$category = factor(coefs$category,
                          levels = c("Outgroup\nRejection", "Political\nPreferences", 
                                     "Ingroup\nAttitudes"))
  
  
  p_eff <- ggplot(coefs, 
              aes(x = dv, y = coef)) +
    geom_hline(yintercept = 0, 
               colour = "red", lty = 2.5) +
    geom_point(aes(x = dv, 
                   y = coef), cex=2) + 
    geom_linerange(aes(x = dv, 
                       ymin = low95,
                       ymax = hi95),
                   linewidth = 1) +
    # geom_linerange(aes(x = dv, 
    #              ymin = low99,
    #              ymax = hi99),
    #            linewidth = .5) + 
    scale_x_discrete(name = "Dependent Variables") +
    scale_y_continuous(name = "ATT") + 
    labs(subtitle = "Effects of New Mosques\n(95% CI)") +
    coord_flip() +
    facet_grid(rows = vars(category),
               scales = "free", space = "free_y") +
    theme_bw()
  
  return(list(p_all = p_all, p_relprop = p_relprop, p_relmean = p_relmean, p_eff = p_eff))
  
}


regs_review <- function(panel_only, shalat = TRUE, how_religious = TRUE, d) {
  
  if (panel_only) {
    
    cols = c("dv", "year", "treat", 
             "female", "age", "edu", "how_religious", "shalat3pt", "work_lastyear",
             "marital", "javaisland", 
             "census10_median_age", "census10_prop_male", "census10_prop_working", "census10_median_edu",
             "census10_prop_muslim", "census10_diversity_lang", 
             "elec14_kec_jokowi_share", "elec19_islamist_share", "elec19_jokowi_share",
             "n_mo_07", 
             "nmprov_mosque", "nmkab_mosque", "nmkec_mosque", "kecid", "pidlink",
             "sharia_pre08")
    
    d = d[complete.cases(d[,cols]),]
    
    ## which pidlink appears twice, meaning they are panel respondents
    vec = d$pidlink[which(duplicated(d$pidlink))]
    
    ## subset to contain only panel respondents
    d = subset(d, subset = pidlink %in% vec)
    
  }
  
  ngrps = length(unique(d$kecid))
  n = length(unique(d$pidlink))
  v = c(ngrps, n)
  names(v) = c("ngrps", "N")
  # print(v)
  
  
  ### basic models
  d$kecid = factor(d$kecid)
  
  if (shalat & how_religious) {
    
    m = lme4::lmer(dv ~ year:treat + year + treat + 
                     female + age + edu + how_religious + shalatlab3pt + work_lastyear + marital + 
                     javaisland + 
                     census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
                     census10_prop_muslim + census10_diversity_lang + 
                     elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
                     ln_n_mo_07 + sharia_pre08 + (1|kecid), 
                   data=d)
  } 
  
  if (!shalat) {
    m = lme4::lmer(dv ~ year:treat + year + treat + 
                     female + age + edu + how_religious + work_lastyear + marital + 
                     javaisland + 
                     census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
                     census10_prop_muslim + census10_diversity_lang + 
                     elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
                     ln_n_mo_07 + sharia_pre08 + (1|kecid), 
                   data=d)
  }
  
  
  if (!how_religious) {
    m = lme4::lmer(dv ~ year:treat + year + treat + 
                     female + age + edu + shalatlab3pt + work_lastyear + marital + 
                     javaisland + 
                     census10_median_age + census10_prop_male + census10_prop_working + census10_median_edu + 
                     census10_prop_muslim + census10_diversity_lang + 
                     elec14_kec_jokowi_share + elec19_islamist_share + elec19_jokowi_share + 
                     ln_n_mo_07 + sharia_pre08 + (1|kecid), 
                   data=d)
  }
  
  m_cl = model_parameters(m, 
                          #effects = "fixed",
                          vcov = "vcovCR",
                          vcov_args = list(type = "CR0", 
                                           cluster = d$kecid))
  
  cols = c("Parameter", "Coefficient", "SE", "t", "p")
  df = data.frame(m_cl)[, cols]
  df$Coefficient = round(df$Coefficient, 4)
  df$SE = round(df$SE, 4)
  df$t = round(df$t, 3)
  df$p = round(df$p, 3)
  
  df = bind_rows(df, data.frame(Parameter = "N Obs", Coefficient = nobs(m)/2, SE = NA, t = NA, p = NA))
  df = bind_rows(df, data.frame(Parameter = "N Districts", Coefficient = ngrps(m), SE = NA, t = NA, p = NA))
  rownames(df) = NULL
  
  df_kab = knitr::kable(df, row.names= TRUE, format = "html", table.attr = "style='width:65%;'")
  
  return(list(m=m, df=df, df_kab = df_kab, data = d))
  
}
